Ion transport through a graphene nanopore 
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Abstract 

Molecular dynamics simulation is utilized to in- 
vestigate the ionic transport of NaCl in solution 
through a graphene nanopore under an applied 
electric field. Results show the formation of con- 
centration polarization layers in the vicinity of the 
graphene sheet. The non-uniformity of the ion dis- 
tribution gives rise to an electric pressure which 
drives vortical motions in the fluid if the electric 
field is sufficiently strong to overcome the influ- 
ence of viscosity and thermal fluctuations. The 
relative importance of hydrodynamic transport and 
thermal fluctuations in determining the pore con- 
ductivity is investigated. A second important ef- 
fect that is observed is the mass transport of wa- 
ter through the nanopore, with an average veloc- 
ity proportional to the applied voltage and inde- 
pendent of the pore diameter. The flux arises as 
a consequence of the asymmetry in the ion distri- 
bution with respect to reflection about the plane 
of the graphene sheet. The accumulation of liq- 
uid molecules in the vicinity of the nanopore due 
to reorientation of the water dipoles by the local 
electric field is seen to result in a local increase in 
the liquid density. Results confirm that the electric 
conductance is proportional to the nanopore diam- 
eter for the parameter regimes that we simulated. 
The occurrence of fluid vortices is found to result 
in an increase in the effective electrical conduc- 
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tance. 

Introduction 

Ionic conduction through nanometer sized chan- 
nels or pores is a common theme in biological sys- 
tems as well as in various manufactured materials 
such as membranes and synthetic nanopores.- - - 
Theoretical as well as experimental aspects of 
the problem have attracted increasing interest in 
the past decade. Molecular dynamics (MD) sim- 
ulation is an effective tool for exploring these 
nano scale phenomena. It has the advantage of be- 
ing able to relate the observables directly to the 
molecular properties of the solid and liquid, once 
an appropriate intermolecular potential is given, 
without the need for too many simplifying as- 
sumptions. Various aspects of the problem, such as 
ionic current rectification, DNA translocation and 
water transport have already been reported in the 
literature.- — 

MD simulations in the context of transloca- 
tion of single stranded and double stranded DNA 
through biological a-Hemolysin and synthetic 
nanopores have been reported.-^*^ It was found 
that the open-pore current increases linearly with 
the applied voltage, and, obeys Ohm's law for volt- 
ages that are no more than of the order of a Volt. 
The distribution of the electric potential around 
the pore, the translocation speed of DNA, the in- 
teractions between the DNA molecule with the 
pore wall, as well as strategies for controlling the 
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transolcation speed have been considered. ' Re- 
cently, Sathe et al.-^ investigated the translocation 
of DNA through a graphene nanopore using MD 
simulations, and suggested that nucleotide pairs 
can be discriminated using graphene nanopores 
under suitable bias conditions. 

There are, however, a number of open problems 
that have not been understood. First, due to the 
geometry of the system, the ions may accumu- 
late and form a concentration polarization layer 
(CPL) in the vinicity of the membrane. The in- 
fluence of the nonuniformity of the ion concentra- 
tion in this charge separated Debye layer on the 
ionic current has not been well understood. Sec- 
ondly, open questions remain on the effect of hy- 
drodynamic flow on the ionic current as well as 
on DNA translocation speeds. Ghosal 15,16 pre- 
sented a simple hydrodynamic calculation for the 
electrophoretic speed of the polymer, modeled lo- 
cally as a long cylindrical object centered on the 
axis of the pore. The calculation yielded analyt- 
ical results in close agreement with experimental 
measurements.— 1 ^ In a related problem where the 
electrophoretic force was measured with the DNA 
immobilized in the pore, numerical 19 as well as 
analytical 20 models based on the hypothesis that 
the hydrodynamic drag was the primary resistive 
force yielded results in close agreement to the ex- 
periments. These results point to the possibility 
that hydrodynamics might play an important role 
in determining DNA translocation speeds. Hydro- 
dynamics also plays an important role in a related 
problem; when a direct current is applied across 
an ion-selective nanoporous membrane or through 
a nanochannel with overlapping Debye layers, it is 
known, that, micro fluid vortices may be observed 
due to hydrodynamic instability, and these vortices 
are capable of enhancing the ionic current in the so 
called "overlimiting regime".— 

The purpose of the present study is to discuss the 
influence of the CPL and fluid convection on ionic 
transport through nanochannels using MD simula- 
tions. The rest of this paper is organized as fol- 
lows: in Sec. 2, we present the physical model and 
describe the numerical approach. The simulation 
results are analysed in detail in Sec. 3. Finally, 
some concluding remarks are made in Sec. 4. 



Method 
System Setup 

The molecular simulations are conducted in a cu- 
bic box with dimensions of L x = 5. 1 12 nm, Ly = 
5.184 nm, L z = 10 nm. The origin of the coordi- 
nates is set in the center of the box. The graphene 
sheet with the size of L x x L y , consists of an ar- 
ray of carbon atoms with a planar hexagonal struc- 
ture (neighboring atoms have a bond lengths of 
0. 142 nm), localized at the mid-plane, z = 0. A 
nanopore, with radius a, is constructed by remov- 
ing the carbon atoms in a central circular patch 
of the graphene sheet: x 2 + y 2 < a 2 . The remain- 
der of the box is then filled with water molecules 
described by the extended simple point charge 
(SPC/E) model. 22 Polarization of water molecules 
is neglected, prior investigations have shown that 
this is reasonable as long as the field strength does 
not exceed about 10 V/nm. 23 Salt (NaCl) is intro- 
duced at a given concentration of 1M by replacing 
the required number of water molecules with Na + 
and Cl~ ions. 

Molecular Dynamics (MD) simulation 

The simulations are performed at a constant tem- 
perature 300K and pressure of 1 bar with the large 
scale MD package GROMACS 4.5. 4.-^ The van 
der Waals (vdW) interaction of the carbon atoms 
is modeled as uncharged Lennard- Jones (LJ) parti- 
cles. The graphene-water interaction is considered 
by a carbon-oxygen LJ potential. This general set 
up and parameter values have been employed in 
previous studies.— - — 

A uniform external electric field is applied in 
directions perpendicular to the graphene sheet (z- 
direction). The contribution of the external elec- 
tric field is described as U E = — ^ 1i r i ' E, where 
the qi and r, denote the charge and location of the 
charged atom i respectively, and E is the strength 
of the external electric field. The LJ interactions 
are truncated at the cut-off distance tq = 1.0 nm 
and the particle mesh Ewald (PME) method 28 with 
a real-space cut-off of 1 nm is utilized to treat 
the long-range electrostatic interactions. Periodic 
boundary conditions are imposed in all directions. 
The time step in all simulations is set to be 2 fs. 
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For the sake of computational efficiency, all of 
the carbon atoms are frozen during the simula- 
tions. Previous investigations have shown that this 
only has a minor influence on the dynamics of the 
adjacent water. Minimization of energy is per- 
formed with the steepest descent method on the 
initial system. Then the system is evolved for 2 ns 
to achieve a state of statistical equilibrium. In all 
cases, statistics are collected during the last 8 ns 
and samples are taken every 0.2 ps. 



Q]for different grid sizes. After determining the av- 



erage velocity and standard error a, it is seen 
that the PDF is well represented by a Gaussian: 
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Furthermore, the PDFs are nearly independent of 
the computational grid size. In the present study, 
we set Az = 0.2 nm. 



Data analysis 

To analyse the statistics of the macroscopic physi- 
cal variables, the spatial location and velocity vec- 
tors of the particles are transformed to cylindrical 
coordinates since averages of variables should be 
close to axisymmetric. Data is presented in the 
r — z plane which is partitioned into gridded cells 
with interval range (Ar,Az ). 

Ionic concentration. The number density 
of Na and CI ions in the cells are calculated by p ; = 
N c /V c , where N c and V c are the number of ions and 
volume of the corresponding cell respectively. 

Water flux. The number of oxygen atoms N w 
through the nanopore in time interval At is cal- 
culated to obtain the water flux. The average wa- 
ter velocity v p in the axial direction through the 
nanopore is given by v p = mN w / AtpoA, in which 
m is the mass of a water molecule, po is water 
density in the bulk and A = Tta 2 is the area of the 
nanopore. 

Flow field. For the cell at (r, z) containing N c 
water molecules, the transient local velocity vector 
can be obtained by 

I N c 

V(r ' z '^ = AtN c ^ [ r oi{tj + At)-roi(tj)] 

where Voiitj) is the location vector of oxygen atom 
i at time tj. At is the time interval between two 
successive frames. We take At = 2 ps in the 
present study. The velocity field v(r,z) of the equi- 
librium state is then calculated by averaging the 
transient velocity over time. To check the indepen- 
dence of physical quantities with respect to grid 
size, the probability density function (PDF) of the 
axial velocity vf of water at (0, 0) is plotted in Fig. 
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Figure 1: Probability distribution function of z- 
component of flow velocity at the nanopore for dif- 
ferent sizes of grid. The velocity distribution can 
be well represented by a Gaussian. 

Ionic current. The ionic current is obtained 
by / = eNi/At, where N, stands for the number of 
Na or CI ions across the graphene nanopore in time 
interval At, e is the proton charge. 

Results and Discussions 

Concentration polarization 

If the pore size approaches zero, clearly a po- 
larized double layer would form adjacent to the 
graphene sheet. For a finite, but nevertheless suf- 
ficiently small pore size, a nonuniform ion distri- 
bution (concentration polarization) is expected ad- 
jacent to the sheet. The profiles of ionic number 
density along the z-axis is plotted in Fig. [2l where 
z = corresponds to the graphene sheet. It is inter- 
esting that the sodium and chlorine ions both form 
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Figure 2: Number density profile of Na and CI ions in z-direction with d=l.5 nm for different field 
strengths, showing the presence of concentration polarization layer. 



concentration polarization layers (CPL) on either 
side of the graphene sheet. However, their distri- 
butions show some asymmetry. Due to size exclu- 
sion, the ions cannot approach the graphene sheet 
closer than about an ionic radius. The combina- 
tion of electrostatic and these steric forces result 
in the appearance of a concentration peak. In the 
case of the chlorine ions, the concentration peak is 
located a distance D c i = 0.375 nm from the sheet. 
This distance is found to be independent of the ap- 
plied field strength. The CPL can also be observed 
for the sodium ions on the other side; however, the 
amplitude of the peak is comparatively weak. For 
weak applied fields (e.g. E = 0.2 V/nm) the CPL 
is not as well defined. A second smaller peak of 
sodium ions may be observed on the other side of 
the membrane. The asymmetry between the two 
kinds of ions may be partly due to the higher mo- 
bility of the chlorine ions, and, partly due to dif- 
ferences in the van der Waals interactions between 
the two kinds of ions with the carbon atoms in 
the graphene sheet. The peak in the ionic concen- 
tration profile is found to increase monotonically 
with increasing electric field. 

In Fig. [31 the number density of CI ions with 
E = 0.5 V/nm and nanopore diameter d = 1.5 nm 
is observed to exhibit a well defined peak. The 
nonuniformity of the ionic distribution will give 



rise to spatial gradients of the electric field strength 
near the membrane, especially in regions close to 
the nanopore. Consequently, variations of the elec- 
tric pressure is to be expected, which might drive 
a flow, as long as these forces are strong enough 
to overcome viscous resistance and are not com- 
pletely masked by the fluctuating Brownian forces 
(thermal fluctuations). 

Flow fields 

The formation of the CPL is expected to generate 
a fluid flow in the vicinity of the interface. This is 
indeed seen in the simulations. The flow stream- 
lines for the nanopore of diameter d = 1.5 nm, 
together with the distribution of water density is 
shown in Fig. |4] with different values of electric 
field strength. When the field is sufficiently large, 
e.g., E = 0.5 and 1.0 V/nm, vortices of a spatial 
scale on the order of nanometers are clearly ob- 
served. When £=0.2 V/nm, the electric pressure 
appears too weak to generate sustained micro vor- 
tices against the influence of thermal fluctuations. 

An interesting physical insight that emerges 
from these simulations is that there appears to 
be a marked asymmetry in the distribution of the 
cations and anions and that of water density. This 
is most likely due to a combination of three fac- 
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Figure 3: Nonuniform distributions of number density of chlorine (a) and sodium (b) ion with d=\.5 nm 
and £=0.5 V/nm. The grey rectangle at z=0 marks the graphene sheet. Ionic gradients lead to electric 
pressure, which drives nanoscale vortices near the pore. 



tors. First, as we have mentioned above, there is 
an approximately 30 percent difference in mobil- 
ity between sodium and chlorine ions. Second, the 
van der Waals interactions between Na and C, CI 
and C differ significantly. The third, which has 
been reported by previous researchers,— - — is that 
the direction of the electric field could have signif- 
icant effects on the interfacial water structures, as 
well as on the wetting behavior of the solid mem- 
brane, due to the polarity of the water molecules. 
A consequence of this asymmetry, is that there is 
a flux of water through the nanopore as may be 
seen in the appearance of the streamlines in Fig. 
IU The average flow velocity v p , describing wa- 
ter crossing the nanopore, is presented in Fig. |5] 
for different electric field strengths and pore diam- 
eters. It shows that generally the average velocity 
at the nanopore is proportional to the applied volt- 
age, and nearly independent of the diameter of the 
nanopore for the parameters we simulated. 

A region of high water density is seen at the cen- 
ter of the nanopore. For E = 1.0 V/nm, the local 
water density can be as high as 2043 kg/m 3 . The 
accumulation of water in the nanopore appears to 
originate from re-orientation of the water dipole 



moments when the external voltage is applied. To 
describe the orientation of water molecules, we 
define an average over the x-y plane of the an- 
gle between the dipole moment vector of the 
water molecule and the z-axis. In the absence 
of an electric field, the orientation of the dipoles 
is completely random corresponding to = 90°. 
Fig. [6] shows that with increasing electric field 
strengths, the dipoles prefer an axial orientation. 
For a given z-location, this effect increases with 
the field strength. The degree of ordering also 
increases as one approaches the pore where the 
dipoles prefer to be parallel to the graphene sur- 
face. This is consistent with previous studies on 
SWCNT where a similar ordering has been ob- 
served. 33 Such ordering allows closer packing of 
the water molecules leading to a rise in density in 
the pore region. 

Voltage current relations 

Recent experimental measurements have produced 
apparently conflicting results on the dependence of 
pore conductance on pore radius. Garaj et al.— re- 
ported that the pore conductance is proportional to 
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Figure 4: Flow fields and water density (in units of kg/m 3 ) distribution of salt solution for a nanopore of 
diameter d=l.5 nm with (a) £=0.2 V/nm; (b) £=0.5 V/nm; (c) £=1.0 V/nm. The rectangles (color grey 
near z=0) marks the graphene sheet. For strong fields, as in panel (c), a well defined vortical flow is seen 
but for weak electric fields the bulk motion of the fluid is difficult to distinguish from the background 
thermal fluctuations. 
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Figure 5: Variation of the mean flow (v p ) at the 
nanopore with the applied electric field for differ- 
ent pore diameters. 

the pore diameter, whereas Schneider et al.— re- 
ported it to be proportional to the pore area. The 
former result would correspond to a hole in an in- 
finitely thin membrane in a homogeneous infinite 
conducting medium. The latter would correspond 



Figure 6: Average angle between the dipole mo- 
ment and the z-axis for a pore diameter of d = 
1.5 nm. Reduction of the angle due to the re- 
orientation of the dipole moment by the field is 
strongest in the pore region where the dipoles fa- 
vor to be parallel to the graphene sheet. 

to a cylindrical conduit of diameter much less than 
the length of the cylinder. 36 Since the former re- 
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E(V/nm) d(nm) 

Figure 7: The current vs applied field (left panel) for several different pore diameters. The pore conduc- 
tance determined from the slope of the initial linear region of the current- voltage characteristic is also 
shown as a function of the pore diameter (right panel). 



suit is to be expected if the thickness of graphene 
can be ignored and the electrolyte solution is as- 
sumed to obey Ohm's law, some researchers 3 have 
speculated that the discrepancy between the ex- 
periment of Schneider et al. and the theoretical 
result for a hole in a thin sheet might originate 
from the fact that their graphene sheet was treated 
with a polymer coating (to reduce non-specific sur- 
face adsorption of DNA). To clarify this, Sathe 
et al.— conducted molecular dynamics simulation 
for a KC1 solution with pore diameters in the range 
of 2 — 7 nm. They found the dependence of re- 
sistance R on pore diameter follows the relation- 
ship R ~ l/d 2 , which is qualitatively in agreement 
with the experiment of Schneider et al.— How- 
ever, quantitatively the resistances they obtained in 
their simulations were three to four times smaller 
than the experimental measurements. They at- 
tributed the discrepancy to a number of factors in- 
cluding inaccuracies in the force field in their sim- 
ulations, unknown charge distribution and uncer- 
tainties about the exact shape of graphene pores. 

The variation of the current with the applied 
electric field was extracted from our simulations 
and is shown in Fig. [7] for a number of differ- 
ent pore diameters, d. A linear regime is ob- 



served at weak fields (E). For higher applied volt- 
ages, it is found that the current shows super-linear 
growth and the nonlinearity is stronger for pores of 
smaller diameter. Hydrodynamic transport due to 
the nanoscale vortices could be responsible for this 
faster than linear increase of the current. A simi- 
lar phenomenon has also been observed when an 
external electric field is applied normal to the sur- 
face of an ion-selective nanoporous membrane im- 
mersed in an electrolyte solution. There, at a crit- 
ical value of the applied field, the quiescent state 
undergoes a hydrodynamic instability resulting in 
the appearance of micro vortices that are respon- 
sible for an "overlimiting-current" in the current- 
voltage relationship. 21 In the nanopore problem 
however, there is no instability. The hydrody- 
namic flow is always present, its strength simply 
increases with the electric field. 

Fig. Hb) shows the electric conductance (G) at 
low fields, which is calculated from the slope of 
the initial linear segment of the I — V curves. It is 
seen to be proportional to the nanopore diameter 
d, which is qualitatively in agreement with the ex- 
periments by Garaj et al. 34 and is as expected from 
the theoretical model where the pore is regarded as 
a hole in an infinitely thin insulating sheet in a uni- 
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form conductor. The resistance obtained for d = 3 
nm is R = 199 MQ, which is comparable to the 
values reported by Schneider et al, 35 though we do 
not observe their G °= d 2 dependence on pore di- 
ameter. The theoretical result G = od for a hole of 
diameter 'J' in an insulating sheet in a conducting 
medium of conductivity 'a' may be used to extract 
an "effective" conductivity from the data shown in 
Fig. |7j5. The conductivity obtained in this way 
is o = 1.879 Sm -1 , which is about a factor of 5 
smaller than the bulk conductivity, o = 9.3 Sm _1 , 
for a 1 M NaCl solution. The discrepancy may 
arise from the fact that there are significant dif- 
ferences between a graphene nanopore immersed 
in an electrolyte and a hole in an insulating sheet 
in an Ohmic medium. As we have seen, one im- 
portant difference lies in the existence of the con- 
centration polarization layer which modifies the 
electric field distribution near the pore as well as 
drives a hydrodynamic flow that affects transport 
properties. The influence of these phenomena on 
pore conductance is still an open problem. A sec- 
ond source of discrepancy may be due to a lim- 
itation of our computation. In a physical experi- 
ment, the nanopore radius is vanishingly small in 
comparison to the reservoir dimensions, so that es- 
sentially all of the electrical resistance arises from 
the pore. This assumption might not be satisfied as 
well in our simulated model, since, due to limita- 
tions of computational capabilities, the area of the 
graphene sheet that we constructed, 5nm x 5nm, 
is quite comparable to the pore diameter of 3 nm. 
These issues will be investigated further in future 
studies. 

Conclusion 

Molecular dynamics simulations were performed 
to study the transport properties of an electrolyte 
through a nanometer sized pore. Due to the par- 
tial blockage of the ionic current by the graphene 
membrane, a concentration polarization layer 
(CPL) develops next to the membrane when an 
external voltage is applied. The CPL is able to 
induce electric pressure in the fluid adjacent to the 
pore. If the applied voltage is large enough to over- 
come the effects of viscosity and Brownian fluctu- 
ations, vortices are generated in the fluid near the 



pore. These nanoscale vortices enhance the ionic 
current through a mechanism similar to an effect 
responsible for the phenomenon of "overlimiting 
current" for perm selective membranes. Owing to 
the polarization of water molecules, the direction 
of the electric field might have significant influ- 
ence on water structure near the graphene surface. 
This effect, together with the differing mobility 
and van der Waals attraction to C exhibited by 
the CI and Na ions, brings about an asymmetric 
distribution of ions, and a net flux of water in the 
direction of the electric field through the pore. 
The average velocity of this hydrodynamic flow is 
found to increase approximately linearly with the 
electric field and appears to be independent of the 
pore diameter. Orientational order created in the 
water dipoles by the electric field enhances closer 
packing of the water dipoles. This is manifested 
in the appearance of regions of high water den- 
sity within the pore. The electrical conductance 
of the system is found to vary linearly with the 
nanopore diameter, as one might expect from the 
classical theoretical relation G = od for a hole 
in an insulating membrane within a conductor. 
However, the effective conductivity, G/d is found 
to be nearly 5 times the bulk conductivity of the 
electrolyte. The linear dependence of G on d is in 
accord with the experiments of Garaj et al. but we 
do not observe the G« d 2 dependence reported by 
Schneider et al. This could be due to differences 
in experimental techniques as a result of which ex- 
perimental conditions in the latter experiment do 
not correspond to the model studied in this paper. 
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